#----- Load posterior samples
load("result_PS.RData")

#----- P.S. plot
pdf("PSPlot.pdf", width=12, height=5)

EAE1 <- result[ ,1:7]
EDE  <- result[ ,8:14]
EAE2 <- result[ ,15:21]

p.EAE1 <- apply(result[ ,c(23,26,29,32,35,38,41)], 2, mean)
p.EDE  <- apply(result[ ,c(22,25,28,31,34,37,40)], 2, mean)
p.EAE2 <- apply(result[ ,c(24,27,30,33,36,39,42)], 2, mean)

plot(apply(EDE, 2, mean), ylim=c(-1, 0.4), xlim=c(1,7.2), col="blue",pch=19, cex=p.EDE*15, xaxt="n", xlab="", ylab="E[ Y(1) - Y(0)]")
axis(1, at=c(1,2,3,4,5,6,7), labels=c("SO2","NOx","CO2","SO2, NOx","SO2, CO2","NOx, CO2","SO2, NOx, CO2"))
points(apply(EAE1, 2, mean),col="red", pch=19, cex=p.EAE1*15)
points(apply(EAE2, 2, mean),col="grey", pch=19, cex=p.EAE2*15)
text(c(1+0.35,2+0.3,3+0.3,4+0.2,5+0.16,6+0.16,7+0.16), apply(EAE1, 2, mean)+c(0,-0.02,0,0,0,0,0), (round(p.EAE1,2)), cex=0.8, col="red")
text(c(1+0.3,2+0.3,3+0.3,4+0.18,5+0.18,6+0.21,7+0.18), apply(EDE, 2, mean)+c(0,0.02,0,0,0,0,0), (round(p.EDE,2)), cex=0.8, col="blue")
text(c(1+0.16,2+0.3,3+0.25,4+0.16,5+0.16,6+0.18,7+0.16), apply(EAE2, 2, mean), (round(p.EAE2,2)), cex=0.8, col="grey")
legend(1.5, 0.4, c("EAE1" , "EDE", "EAE2"), pch=c(19,19,19),col=c("red", "blue","grey"))
dev.off()

